ggplot(plot, aes(x=time))+
geom_ribbon(aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u), alpha = 0.5, fill='lightgrey')+
geom_ribbon(data=narrow,aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u), alpha = 0.5, fill='grey')+
geom_line(data=subset(plot, miss_unit!='All'),aes(y=estimated_ATT, color= miss_unit),linewidth=0.5,alpha=0.5)+
geom_line(data=subset(plot, miss_unit=='All'),aes(y=estimated_ATT, color='ATT with full pool'),linewidth=1,alpha=1)+
geom_vline(xintercept=0)+theme_minimal() + theme(panel.background = element_rect(fill = "white"))+
labs(x="Months since treatment", y='ATT with widest and narrowest 95% CIs')+
scale_color_discrete(name='Key',
breaks=c('ATT with full pool', unique(plot$miss_unit)))
ggsave('output/loo_effect.pdf', width = 10, height=5)
ggsave('output/loo_effect.png', width = 10, height=5)
######## SETUP ########
rm(list = ls()) # clear workspace
need <- c("Synth", "reshape2","devtools", "synthdid",'tidyverse',
'bpCausal', "fixest","fwildclusterboot","modelsummary",
'kableExtra') # list packages needed
have <- need %in% rownames(installed.packages()) # checks packages you have
if(any(!have)) install.packages(need[!have]) # install missing packages
#devtools::install_github("synth-inference/synthdid") #manually install sdid if the codes above fail to do so automatically
invisible(lapply(need, library, character.only=T)) # load needed packages
setwd(gsub('/code','',dirname(rstudioapi::getSourceEditorContext()$path))) #set wd to root directory
list.files()
##################DM-LFM####################################
load('stored/dmlfm_import.rdata')
source('custom_functions/DM-LFM.r')
result.use <- effSummary(out) #extract effect information
class(result.use) <- 'DMLFM'
#plot.DMLFM(out, ci.level=0.90)
result.use90 <- effSummary(out, ci.level = 0.90) #for 90% CI
class(result.use90) <- 'DMLFM'
#trend plot
ggplot(result.use$est.eff, aes(x=time))+
geom_ribbon(aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u, fill = "95% CI"), alpha = 0.5)+
geom_ribbon(data=result.use90$est.eff,aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u, fill = "90% CI"), alpha = 0.5)+
geom_line(aes(y=observed, color='Observed'),linewidth=0.75,alpha=0.5)+
geom_line(aes(y=estimated_counterfactual, color='Counterfactual'),linewidth=0.75,alpha=0.5)+
geom_vline(xintercept=0)+theme_minimal() + theme(panel.background = element_rect(fill = "white"))+
scale_color_manual(name= '',
breaks=c('Observed','Counterfactual'),
values=c('Observed'="#00BFC4", "Counterfactual"="#F8766D"))+
labs(x="Months since treatment", y='Log(Filipino imports from China)')+
scale_fill_manual(name=c("", ''),
breaks=c('95% CI', '90% CI'),
values=c('lightgrey', 'grey'))
ggsave('output/dmlfm_trend_import.pdf', width = 10, height=5)
ggsave('output/dmlfm_trend_import.png', width = 10, height=5)
#effect plot
ggplot(result.use$est.eff, aes(x=time, y=estimated_ATT))+
geom_ribbon(aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u, fill = "95% CI"), alpha = 0.5)+
geom_ribbon(data=result.use90$est.eff,aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u, fill = "90% CI"), alpha = 0.5)+
geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +
theme_bw()+labs(x="Months since treatment",y="ATT with 90% and 95% CIs") + geom_line(linewidth=0.75) +
scale_fill_manual(name=c("", ''),
breaks=c('95% CI', '90% CI'),
values=c('lightgrey', 'grey'))         #effect plot
ggsave('output/dmlfm_effect_import.pdf', width = 10, height=5)
ggsave('output/dmlfm_effect_import.png', width = 10, height=5)
result.use$obs <- nrow(data)
result.use$tr <- length(unique(data$Country[data$treated==1]))
result.use$ct <- length(unique(data$Country))-result.use$tr
class(result.use) <- 'DMLFM'
result.use90$obs <- nrow(data)
result.use90$tr <- length(unique(data$Country[data$treated==1]))
result.use90$ct <- length(unique(data$Country))-result.use90$tr
class(result.use90) <- 'DMLFM'
result.use_import <- result.use
result.use90_import <- result.use90
load('stored/dmlfm.rdata')
result.use <- effSummary(out) #extract effect information
class(result.use) <- 'DMLFM'
result.use90 <- effSummary(out, ci.level = 0.90) #for 90% CI
result.use$obs <- nrow(data)
result.use$tr <- length(unique(data$Country[data$treated==1]))
result.use$ct <- length(unique(data$Country))-result.use$tr
class(result.use) <- 'DMLFM'
result.use90$obs <- nrow(data)
result.use90$tr <- length(unique(data$Country[data$treated==1]))
result.use90$ct <- length(unique(data$Country))-result.use90$tr
class(result.use90) <- 'DMLFM'
tidy.DMLFM <- function(x, ...) {
ret <- data.frame(
term      = 'ATT',
estimate  = x$est.avg[[1]],
conf.low  = x$est.avg[[2]],
conf.high = x$est.avg[[3]],
p.value   = x$est.avg_p)
ret
}
glance.DMLFM <- function(x, ...) {
ret <- data.frame(
nobs  = x$obs,
tr  = x$tr,
ct  = x$ct)
ret
}
gof_map <- list(list("raw" = "nobs", "clean" = "Observations", "fmt" = 0),
list("raw" = "tr", "clean" = "Treated Units", "fmt" = 0),
list("raw" = "ct", "clean" = "Control Units", "fmt" = 0))
modelsummary(list('95% CI'=result.use, '90% CI'=result.use90,
'95% CI'=result.use_import, '90% CI'=result.use90_import),
estimate = "{estimate}", #output='latex', #remove hashtag for latex output
statistic = "[{conf.low}, {conf.high}]",
gof_map = gof_map, escape=T,
notes=c('Equal-tailed Credible Intervals in square brackets.'),
title = 'Main Analysis Result table for the DM-LFM.')  %>%
add_header_above(c(' '=1, 'Exports to China' = 2, 'Imports from China' = 2), bold=T)
load('stored/dmlfm.rdata')
result.use <- effSummary(out) #extract effect information
class(result.use) <- 'DMLFM'
result.use90 <- effSummary(out, ci.level = 0.90) #for 90% CI
result.use$obs <- nrow(data)
result.use$tr <- length(unique(data$Country[data$treated==1]))
result.use$ct <- length(unique(data$Country))-result.use$tr
class(result.use) <- 'DMLFM'
result.use90$obs <- nrow(data)
result.use90$tr <- length(unique(data$Country[data$treated==1]))
result.use90$ct <- length(unique(data$Country))-result.use90$tr
class(result.use90) <- 'DMLFM'
modelsummary(list('95% CI'=result.use, '90% CI'=result.use90,
'95% CI'=result.use_import, '90% CI'=result.use90_import),
estimate = "{estimate}", #output='latex', #remove hashtag for latex output
statistic = "[{conf.low}, {conf.high}]",
gof_map = gof_map, escape=T,
notes=c('Equal-tailed Credible Intervals in square brackets.'),
title = 'Main Analysis Result table for the DM-LFM.')  %>%
add_header_above(c(' '=1, 'Exports to China' = 2, 'Imports from China' = 2), bold=T)
##################DM-LFM####################################
data <- read.csv("input/export_to_China.csv")
#logimport
data$logimport <- log(data$C.import)
data$treated <- ifelse(data$Unit == 1 & data$M.unit >= 0, 1,0 ) #data cleaning
set.seed(666)
out = bpCausal(data = data,
Yname = "logimport",
Dname = "treated",
Xname = c(),
Zname = c(),
Aname = c(),
index = c("Country", "M.unit"),
re = "both",
ar1 = TRUE,
r = 10,
niter = 15000,
burn = 5000,
xlasso = 1,
zlasso = 1,
alasso = 1,
flasso = 1,
a1 = 0.001,
a2 = 0.001,
b1 = 0.001,
b2 = 0.001,
c1 = 0.001,
c2 = 0.001,
p1 = 0.001,
p2 = 0.001)  #DM-LFM Warning: a bit time consuming
save(out, data, file='stored/dmlfm.rdata')  #save model (you can load it next time)
load('stored/dmlfm.rdata')
result.use <- effSummary(out) #extract effect information
class(result.use) <- 'DMLFM'
result.use90 <- effSummary(out, ci.level = 0.90) #for 90% CI
result.use$obs <- nrow(data)
result.use$tr <- length(unique(data$Country[data$treated==1]))
result.use$ct <- length(unique(data$Country))-result.use$tr
class(result.use) <- 'DMLFM'
result.use90$obs <- nrow(data)
result.use90$tr <- length(unique(data$Country[data$treated==1]))
result.use90$ct <- length(unique(data$Country))-result.use90$tr
class(result.use90) <- 'DMLFM'
modelsummary(list('95% CI'=result.use, '90% CI'=result.use90,
'95% CI'=result.use_import, '90% CI'=result.use90_import),
estimate = "{estimate}", #output='latex', #remove hashtag for latex output
statistic = "[{conf.low}, {conf.high}]",
gof_map = gof_map, escape=T,
notes=c('Equal-tailed Credible Intervals in square brackets.'),
title = 'Main Analysis Result table for the DM-LFM.')  %>%
add_header_above(c(' '=1, 'Exports to China' = 2, 'Imports from China' = 2), bold=T)
modelsummary(list('95% CI'=result.use, '90% CI'=result.use90,
'95% CI'=result.use_import, '90% CI'=result.use90_import),
estimate = "{estimate}", #output='latex', #remove hashtag for latex output
statistic = "[{conf.low}, {conf.high}]",
gof_map = gof_map, escape=T,
notes=c('Equal-tailed Credible Intervals in square brackets.'),
title = 'Main Analysis Result table for the DM-LFM.')  %>%
add_header_above(c(' '=1, 'Exports to China' = 2, 'Imports from China' = 2), bold=T)
#Short term result table
class(result.use) <- 'DMLFMtime'
class(result.use90) <- 'DMLFMtime'
class(result.use_import) <- 'DMLFMtime'
class(result.use90_import) <- 'DMLFMtime'
tidy.DMLFMtime <- function(x, ...) {
ret <- data.frame(
term      = x$est.eff$time,
estimate  = x$est.eff$estimated_ATT,
conf.low  = x$est.eff$estimated_ATT_ci_l,
conf.high = x$est.eff$estimated_ATT_ci_u)
ret
}
glance.DMLFMtime <- function(x, ...) {
ret <- data.frame(
nobs  = x$obs,
tr  = x$tr,
ct  = x$ct)
ret
}
gof_map <- list(list("raw" = "nobs", "clean" = "Observations", "fmt" = 0),
list("raw" = "tr", "clean" = "Treated Units", "fmt" = 0),
list("raw" = "ct", "clean" = "Control Units", "fmt" = 0))
coef_map <- list('0' = '0',
'1' = '1',
'2' = '2',
'3' = '3',
'4' = '4',
'5' = '5',
'6' = '6')
modelsummary(list('95% CI'=result.use, '90% CI'=result.use90,
'95% CI'=result.use_import, '90% CI'=result.use90_import),
estimate = "{estimate}", #output='latex',
statistic = "[{conf.low}, {conf.high}]",
gof_map = gof_map, coef_map=coef_map, escape=T,
notes=c('Equal-tailed Credible Intervals in square brackets.'),
title = 'Result table for the ATT over time.') %>%
add_header_above(c(' '=1, 'Exports to China' = 2, 'Imports from China' = 2), bold=T)
#Placebo test (DMLFM)
load('stored/dmlfm_import_placebo.rdata')
##################time placebo###############################
data <- read.csv("input/import_from_China.csv")
data_placebo <- subset(data, M.unit<0)
data_placebo$logexport <- log(data_placebo$C.export)
data_placebo$treated <- ifelse(data_placebo$Unit == 1 & data_placebo$M.unit >= -24, 1,0 ) #data cleaning for synthdid
set.seed(666)
out = bpCausal(data = data_placebo,
Yname = "logexport",
Dname = "treated",
Xname = c(),
Zname = c(),
Aname = c(),
index = c("Country", "M.unit"),
re = "both",
ar1 = TRUE,
r = 10,
niter = 15000,
burn = 5000,
xlasso = 1,
zlasso = 1,
alasso = 1,
flasso = 1,
a1 = 0.001,
a2 = 0.001,
b1 = 0.001,
b2 = 0.001,
c1 = 0.001,
c2 = 0.001,
p1 = 0.001,
p2 = 0.001)  #DM-LFM Warning: a bit time consuming
save(out, data_placebo, file='stored/dmlfm_import_placebo.rdata')  #save model (you can load it next time)
#Placebo test (DMLFM)
load('stored/dmlfm_import_placebo.rdata')
result.use <- effSummary(out) #extract effect information
class(result.use) <- 'DMLFM'
result.use90 <- effSummary(out, ci.level = 0.90) #for 90% CI
#trend plot
ggplot(result.use$est.eff, aes(x=time-25))+
geom_ribbon(aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u, fill = "95% CI"), alpha = 0.5)+
geom_ribbon(data=result.use90$est.eff,aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u, fill = "90% CI"), alpha = 0.5)+
geom_line(aes(y=observed, color='Observed'),linewidth=0.75,alpha=0.5)+
geom_line(aes(y=estimated_counterfactual, color='Counterfactual'),linewidth=0.75,alpha=0.5)+
geom_vline(xintercept=-25)+theme_minimal() + theme(panel.background = element_rect(fill = "white"))+
scale_color_manual(name= '',
breaks=c('Observed','Counterfactual'),
values=c('Observed'="#00BFC4", "Counterfactual"="#F8766D"))+
labs(x="Months since actual treatment", y='Log(Filipino imports from China)')+
scale_fill_manual(name=c("", ''),
breaks=c('95% CI', '90% CI'),
values=c('lightgrey', 'grey'))
ggsave('output/placebo_trend_import.pdf', width = 10, height=5)
ggsave('output/placebo_trend_import.png', width = 10, height=5)
#effect plot (90% and 95% CIs)
ggplot(result.use$est.eff, aes(x=time-25, y=estimated_ATT))+
geom_ribbon(aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u), fill = "lightgrey", alpha = 0.5)+
geom_ribbon(data=result.use90$est.eff,aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u), fill = "grey", alpha = 0.5)+
geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
geom_vline(xintercept = -25, linetype = "dashed", color = "grey") +
geom_line(linewidth=0.75)+
theme_bw()+labs(x="Months since treatment",y="ATT with 90% and 95% CIs")
result.use$tr <- length(unique(data$Country[data$treated==1]))
result.use$ct <- length(unique(data$Country))-result.use$tr
class(result.use) <- 'DMLFM'
result.use90$obs <- nrow(data)
result.use90$tr <- length(unique(data$Country[data$treated==1]))
result.use90$ct <- length(unique(data$Country))-result.use90$tr
class(result.use90) <- 'DMLFM'
result.use_import <- result.use
result.use90_import <- result.use90
load('stored/dmlfm_placebo.rdata')
result.use <- effSummary(out) #extract effect information
class(result.use) <- 'DMLFM'
result.use90 <- effSummary(out, ci.level = 0.90) #for 90% CI
result.use$obs <- nrow(data)
result.use$tr <- length(unique(data$Country[data$treated==1]))
result.use$ct <- length(unique(data$Country))-result.use$tr
class(result.use) <- 'DMLFM'
result.use90$obs <- nrow(data)
result.use90$tr <- length(unique(data$Country[data$treated==1]))
result.use90$ct <- length(unique(data$Country))-result.use90$tr
class(result.use90) <- 'DMLFM'
modelsummary(list('95% CI'=result.use, '90% CI'=result.use90,
'95% CI'=result.use_import, '90% CI'=result.use90_import),
estimate = "{estimate}", #output='latex',
statistic = "[{conf.low}, {conf.high}]",
gof_map = gof_map, escape=T,
notes=c('Equal-tailed Credible Intervals in square brackets.'),
title = 'Placebo result table for the DM-LFM.') %>%
add_header_above(c(' '=1, 'Exports to China' = 2, 'Imports from China' = 2), bold=T)
View(result.use)
View(result.use_import)
View(result.use90)
View(result.use90_import)
result.use$obs <- nrow(data)
result.use$tr <- length(unique(data$Country[data$treated==1]))
result.use$ct <- length(unique(data$Country))-result.use$tr
class(result.use) <- 'DMLFM'
result.use90$obs <- nrow(data)
result.use90$tr <- length(unique(data$Country[data$treated==1]))
result.use90$ct <- length(unique(data$Country))-result.use90$tr
class(result.use90) <- 'DMLFM'
result.use_import <- result.use
result.use90_import <- result.use90
load('stored/dmlfm_placebo.rdata')
result.use <- effSummary(out) #extract effect information
class(result.use) <- 'DMLFM'
result.use90 <- effSummary(out, ci.level = 0.90) #for 90% CI
result.use$obs <- nrow(data)
result.use$tr <- length(unique(data$Country[data$treated==1]))
result.use$ct <- length(unique(data$Country))-result.use$tr
class(result.use) <- 'DMLFM'
result.use90$obs <- nrow(data)
result.use90$tr <- length(unique(data$Country[data$treated==1]))
result.use90$ct <- length(unique(data$Country))-result.use90$tr
class(result.use90) <- 'DMLFM'
modelsummary(list('95% CI'=result.use, '90% CI'=result.use90,
'95% CI'=result.use_import, '90% CI'=result.use90_import),
estimate = "{estimate}", #output='latex',
statistic = "[{conf.low}, {conf.high}]",
gof_map = gof_map, escape=T,
notes=c('Equal-tailed Credible Intervals in square brackets.'),
title = 'Placebo result table for the DM-LFM.') %>%
add_header_above(c(' '=1, 'Exports to China' = 2, 'Imports from China' = 2), bold=T)
load('stored/dmlfm_import_placebo.rdata')
result.use <- effSummary(out) #extract effect information
class(result.use) <- 'DMLFM'
result.use90 <- effSummary(out, ci.level = 0.90) #for 90% CI
result.use$obs <- nrow(data)
result.use$tr <- length(unique(data$Country[data$treated==1]))
result.use$ct <- length(unique(data$Country))-result.use$tr
class(result.use) <- 'DMLFM'
result.use90$obs <- nrow(data)
result.use90$tr <- length(unique(data$Country[data$treated==1]))
result.use90$ct <- length(unique(data$Country))-result.use90$tr
class(result.use90) <- 'DMLFM'
result.use_import <- result.use
result.use90_import <- result.use90
load('stored/dmlfm_placebo.rdata')
result.use <- effSummary(out) #extract effect information
class(result.use) <- 'DMLFM'
result.use90 <- effSummary(out, ci.level = 0.90) #for 90% CI
result.use$obs <- nrow(data)
result.use$tr <- length(unique(data$Country[data$treated==1]))
result.use$ct <- length(unique(data$Country))-result.use$tr
class(result.use) <- 'DMLFM'
result.use90$obs <- nrow(data)
result.use90$tr <- length(unique(data$Country[data$treated==1]))
result.use90$ct <- length(unique(data$Country))-result.use90$tr
class(result.use90) <- 'DMLFM'
modelsummary(list('95% CI'=result.use, '90% CI'=result.use90,
'95% CI'=result.use_import, '90% CI'=result.use90_import),
estimate = "{estimate}", #output='latex', #remove hashtag for latex output
statistic = "[{conf.low}, {conf.high}]",
gof_map = gof_map, escape=T,
notes=c('Equal-tailed Credible Intervals in square brackets.'),
title = 'Placebo result table for the DM-LFM.') %>%
add_header_above(c(' '=1, 'Exports to China' = 2, 'Imports from China' = 2), bold=T)
#leave-one-out test
load('stored/leave-one-out_import.rdata')
narrow <- plot %>%
group_by(time) %>%
summarize(counterfactual_ci_l = max(counterfactual_ci_l), counterfactual_ci_u = min(counterfactual_ci_u))
#trend
ggplot(plot, aes(x=time))+
geom_ribbon(aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u), alpha = 0.5, fill='lightgrey')+
geom_ribbon(data=narrow,aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u), alpha = 0.5, fill='grey')+
geom_line(data=subset(plot, miss_unit!='All'),aes(y=estimated_counterfactual, color= miss_unit),linewidth=0.5,alpha=0.5)+
geom_line(data=subset(plot, miss_unit=='All'),aes(y=estimated_counterfactual, color='Counterfactual (full pool)'),linewidth=1,alpha=1)+
geom_line(aes(y=observed, color='Observed'),linewidth=1,alpha=1)+
geom_vline(xintercept=0)+theme_minimal() + theme(panel.background = element_rect(fill = "white"))+
labs(x="Months since treatment", y='Log(Filipino imports from China)')+
scale_color_discrete(name='Key',
breaks=c('Observed', 'Counterfactual (full pool)', unique(plot$miss_unit)))
ggsave('output/loo_trend_import.pdf', width = 10, height=5)
ggsave('output/loo_trend_import.png', width = 10, height=5)
#effect
narrow <- plot %>%
group_by(time) %>%
summarize(estimated_ATT_ci_l = max(estimated_ATT_ci_l), estimated_ATT_ci_u = min(estimated_ATT_ci_u))
ggplot(plot, aes(x=time))+
geom_ribbon(aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u), alpha = 0.5, fill='lightgrey')+
geom_ribbon(data=narrow,aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u), alpha = 0.5, fill='grey')+
geom_line(data=subset(plot, miss_unit!='All'),aes(y=estimated_ATT, color= miss_unit),linewidth=0.5,alpha=0.5)+
geom_line(data=subset(plot, miss_unit=='All'),aes(y=estimated_ATT, color='ATT with full pool'),linewidth=1,alpha=1)+
geom_vline(xintercept=0)+theme_minimal() + theme(panel.background = element_rect(fill = "white"))+
labs(x="Months since treatment", y='ATT with widest and narrowest 95% CIs')+
scale_color_discrete(name='Key',
breaks=c('ATT with full pool', unique(plot$miss_unit)))
ggsave('output/loo_effect_import.pdf', width = 10, height=5)
ggsave('output/loo_effect_import.png', width = 10, height=5)
#SC result table
load('stored/sc_import.rdata')
sc_estimate_import <- sc_estimate
load('stored/sc.rdata')
tidy.sc_estimate <- function(x, ...) {
ret <- data.frame(
term      = 'ATT',
estimate  = x$estimate,
std.err   = x$se,
p.value   = x$p,
conf.low  = x$estimate+x$se*qnorm(0.025),
conf.high  = x$estimate+x$se*qnorm(0.975))
ret
}
glance.sc_estimate <- function(x, ...) {
ret <- data.frame(
nobs  = 891,
tr  = 1,
ct  = 8)
ret
}
gof_map <- list(list("raw" = "nobs", "clean" = "Observations", "fmt" = 0),
list("raw" = "tr", "clean" = "Treated Units", "fmt" = 0),
list("raw" = "ct", "clean" = "Control Units", "fmt" = 0))
modelsummary(list('Exports to China' = sc_estimate, 'Imports from China' = sc_estimate_import),
estimate = "{estimate}{stars} ({p.value})",
statistic = "[{conf.low}, {conf.high}]", #output='latex',
gof_map = gof_map, escape=T,
notes=c('p-value in parentheses; 95% CI in square brackets.'),
title = 'Result table for the ATT estimated by SCM.')
rm(list = ls()) # clear workspace
need <- c('coda','bpCausal','tidyverse','kableExtra') # list packages needed
have <- need %in% rownames(installed.packages()) # checks packages you have
if(any(!have)) install.packages(need[!have]) # install missing packages
#devtools::install_github("synth-inference/synthdid") #manually install sdid if the codes above fail to do so automatically
invisible(lapply(need, library, character.only=T)) # load needed packages
setwd(gsub('/code','',dirname(rstudioapi::getSourceEditorContext()$path))) #set wd to root directory
list.files()
source('custom_functions/DM-LFM.r')
######## SUMMARY STATS########
load('stored/dmlfm_import.rdata')
result.use_import <- effSummary(out) #extract effect information
load('stored/dmlfm.rdata')
result.use <- effSummary(out) #extract effect information
mc_att <- mcmc(result.use$eff_avg_i)
mc_att_import <- mcmc(result.use_import$eff_avg_i)
sum <- rbind(t(summary(mc_att)$statistics), t(summary(mc_att_import)$statistics))
sum <- cbind('Trade'=c('Exports','Imports'),as.data.frame(sum))
kbl(sum, position='h', centering=T, digits=6,
caption = 'Summary table for MCMC.')
######## QUANTILE STATS########
quant <- rbind(t(summary(mc_att)$quantiles), t(summary(mc_att_import)$quantiles))
quant <- cbind('Trade'=c('Exports','Imports'),as.data.frame(quant))
kbl(quant, position='h', centering=T, digits=6,
caption = 'Quantiles.')
######## TRACE PLOTS #########
png('output/MCMCtraceplot.png', units='in', width=8, height=4, res=300)
traceplot(mc_att)
dev.off()
pdf('output/MCMCtraceplot.pdf', width=8, height=4)
traceplot(mc_att)
dev.off()
png('output/MCMCtraceplot_import.png', units='in', width=8, height=4, res=300)
traceplot(mc_att_import)
dev.off()
pdf('output/MCMCtraceplot_import.pdf', width=8, height=4)
traceplot(mc_att_import)
dev.off()
######## DENSITY PLOTS #########
png('output/MCMCdensplot.png', units='in', width=8, height=4, res=300)
densplot(mc_att)
dev.off()
pdf('output/MCMCdensplot.pdf', width=8, height=4)
densplot(mc_att)
dev.off()
png('output/MCMCdensplot_import.png', units='in', width=8, height=4, res=300)
densplot(mc_att_import)
dev.off()
pdf('output/MCMCdensplot_import.pdf', width=8, height=4)
densplot(mc_att_import)
dev.off()
######## GEWEKE PLOTS #########
geweke.diag(mc_att)
pnorm(geweke.diag(mc_att)$z, lower.tail=T)*2
geweke.diag(mc_att_import)
pnorm(geweke.diag(mc_att_import)$z, lower.tail=F)*2
png('output/MCMCgeweke.png', units='in', width=8, height=4, res=300)
geweke.plot(mc_att, nbins=21)
dev.off()
pdf('output/MCMCgeweke.pdf', width=8, height=4)
geweke.plot(mc_att, nbins=21)
dev.off()
png('output/MCMCgeweke_import.png', units='in', width=8, height=4, res=300)
geweke.plot(mc_att_import, nbins=21)
dev.off()
pdf('output/MCMCgeweke_import.pdf', width=8, height=4)
